Ill-conditioned System Of Linear Equations Research Articles

Full operator preconditioning and the accuracy of solving linear systems

Abstract Unless special conditions apply, the attempt to solve ill-conditioned systems of linear equations with standard numerical methods leads to uncontrollably high numerical error and often slow convergence of an iterative solver. In many cases, such systems arise from the discretization of operator equations with a large number of discrete variables and the ill-conditioning is tackled by means of preconditioning. A key observation in this paper is the sometimes overlooked fact that while traditional preconditioning effectively accelerates convergence of iterative methods, it generally does not improve the accuracy of the solution. Nonetheless, it is sometimes possible to overcome this barrier: accuracy can be improved significantly if the equation is transformed before discretization, a process we refer to as full operator preconditioning (FOP). We highlight that this principle is already used in various areas, including second kind integral equations and Olver–Townsend’s spectral method. We formulate a sufficient condition under which high accuracy can be obtained by FOP. We illustrate this for a fourth order differential equation which is discretized using finite elements.

A multigrid preconditioner for spatially adaptive high-order meshless method on fluid–solid interaction problems

We present a monolithic geometric multigrid preconditioner for solving fluid–solid interaction problems in Stokes limit. The problems are discretized by a spatially adaptive high-order meshless method, the generalized moving least squares (GMLS) with adaptive h-refinement. For solving fluid–solid interaction problems, we need to deal with a tightly coupled system consisting of the flow field and solid bodies, resulting in a linear system of equations with a block structure. In Stokes limit, solid kinematics can be dominated by the singularities governing the lubrication effects. Resolving those singularities with adaptive h-refinement can lead to an ill-conditioned linear system of equations. The key ingredients of the multigrid preconditioner include the interpolation and restriction operators and the smoothers. For constructing the interpolation and restriction operators, we utilize the geometric information of hierarchical sets of GMLS nodes generated in adaptive h-refinement. We build decoupled smoothers through physics-based splitting and then combine them via a multiplicative overlapping Schwarz approach. Through numerical examples with the inclusion of different numbers and shapes of solid bodies, we demonstrate the performance and assess the scalability of the designed preconditioner. As the total degrees of freedom and the number of solid bodies increase, the proposed monolithic geometric multigrid preconditioner can ensure convergence and good scalability when using the Krylov iterative method for solving the linear systems of equations generated from the spatially adaptive GMLS discretization. More specifically, for a fixed number of solid bodies, as the discretization resolution is incrementally refined, the number of iterations of the linear solver can be maintained at the same level, indicating nearly linear scalability of our preconditioner with respect to the total degrees of freedom. When the number of solid bodies Ns increases, the number of iterations is nearly proportional to Ns, implying the sublinear optimality with respect to the number of solid bodies.

Diffusion source function identification in the process of gas seepage in a coal seam

Through piecewise linearization, the control equation of gas migration in coal is simplified to a parabolic differential equation with time-dependent source terms, and under the corresponding initial conditions, boundary and additional conditions are used to obtain a definite solution to the unknown source function identification of the parabolic differential equation. A triaxial gas migration experimental device to measure the axial gas pressure of a specimen is independently developed, and a gas migration test of coal samples is performed. With the use of the fundamental equation solution, by substituting the definite solution conditions obtained in the experiment, the above equation is discretized into a system of linear equations, the Tikhonov regularization method and generalized cross-validation (GCV) method are employed to solve the obtained ill-conditioned linear system of equations, and the diffusion source functions in the seepage equation are identified. The results suggest that the unknown function identification method for the differential equation, which avoids measurement and analysis of the pore and fissure structures of coal samples, accurately and directly obtains the change trend of desorption–diffusion sources in the process of gas seepage and provides a new idea for the study of the coal seam gas migration process.

Maximum Entropy Technique and Regularization Functional for Determining the Pharmacokinetic Parameters in DCE-MRI

This paper aims to solve the arterial input function (AIF) determination in dynamic contrast-enhanced MRI (DCE-MRI), an important linear ill-posed inverse problem, using the maximum entropy technique (MET) and regularization functionals. In addition, estimating the pharmacokinetic parameters from a DCE-MR image investigations is an urgent need to obtain the precise information about the AIF–the concentration of the contrast agent on the left ventricular blood pool measured over time. For this reason, the main idea is to show how to find a unique solution of linear system of equations generally in the form of y=Ax+b, named an ill-conditioned linear system of equations after discretization of the integral equations, which appear in different tomographic image restoration and reconstruction issues. Here, a new algorithm is described to estimate an appropriate probability distribution function for AIF according to the MET and regularization functionals for the contrast agent concentration when applying Bayesian estimation approach to estimate two different pharmacokinetic parameters. Moreover, by using the proposed approach when analyzing simulated and real datasets of the breast tumors according to pharmacokinetic factors, it indicates that using Bayesian inference—that infer the uncertainties of the computed solutions, and specific knowledge of the noise and errors—combined with the regularization functional of the maximum entropy problem, improved the convergence behavior and led to more consistent morphological and functional statistics and results. Finally, in comparison to the proposed exponential distribution based on MET and Newton’s method, or Weibull distribution via the MET and teaching–learning-based optimization (MET/TLBO) in the previous studies, the family of Gamma and Erlang distributions estimated by the new algorithm are more appropriate and robust AIFs.

A stabilized GMRES method for singular and severely ill-conditioned systems of linear equations

Consider using the right-preconditioned GMRES (AB-GMRES) for obtaining the minimum-norm solution of inconsistent underdetermined systems of linear equations. Morikuni (Ph.D. thesis, 2013) showed that for some inconsistent and ill-conditioned problems, the iterates may diverge. This is mainly because the Hessenberg matrix in the GMRES method becomes very ill-conditioned so that the backward substitution of the resulting triangular system becomes numerically unstable. We propose a stabilized GMRES based on solving the normal equations corresponding to the above triangular system using the standard Cholesky decomposition. This has the effect of shifting upwards the tiny singular values of the Hessenberg matrix which lead to an inaccurate solution. We analyze why the method works. Numerical experiments show that the proposed method is robust and efficient, not only for applying AB-GMRES to underdetermined systems, but also for applying GMRES to severely ill-conditioned range-symmetric systems of linear equations.

Measurement of linear accelerator spectra, reconstructed from percentage depth dose curves by neural networks

Purpose Linear accelerator (linac) spectra, used to improve the accuracy of dose calculation and to produce a complete description of beam quality, among other aspects, are relevant in radiotherapy and linear accelerator physics.Methods In this work we apply neural networks in solving an ill-conditioned system of linear equations, to indirectly measure the linear accelerator spectra via the percentage depth dose curves. The photon beam spectra are related to radiation doses through a Fredholm integral equation. To address our problem we measured the percentage depth dose curve in water and solved a discretized Fredholm equation using artificial neural network. After analysing the typology of our physical problem we selected a MultiLayer Perceptron Neural Network and designed the most suitable neural network architecture.Results The reconstructed spectra were compared with spectra from three linacs, two of them obtained by us through simulations and the third produced by another author, achieving a concordance between 92 % and 96 %.Conclusions Therefore, the spectrum of any accelerator can be quickly and easily reconstructed from measured percent depth dose curves, applying a trained artificial neural network.

A recycling preconditioning method with auxiliary tip subspace for elastic crack propagation simulation using XFEM

In this paper, we propose a recycling preconditioning method with auxiliary tip subspace for solving a sequence of highly ill-conditioned linear systems of equations of different sizes arising from elastic crack propagation problems discretized by the extended finite element method. To construct a Schwarz type preconditioner, the finite element mesh is decomposed into crack tip subdomains, which contain all the degrees of freedom (DOFs) of the branch enrichment functions, and regular subdomains, which contain the standard DOFs and the DOFs of the Heaviside and the Junction enrichment functions. As cracks propagate these subdomains are modified accordingly, and the subdomain matrices are constructed as the restriction of the global matrix to the subdomains. In the overlapping Schwarz preconditioners, the crack tip subproblems are solved exactly and the regular subproblems are solved by some inexact solvers, such as ILU. We consider problems with and without crack intersections and develop a simple scheme to update, instead of re-computing, the subdomain problems as cracks propagate, in which only crack tip subdomains are updated around the new crack tips and all the regular subdomains remain unchanged. Therefore, no extra search is required, and the sizes of crack tip subproblems do not increase as cracks propagate, which greatly saves the computational cost. Moreover, starting from the second system, the Krylov subspace method uses a nontrivial initial guess constructed using the solution of the previous system with a modification around the new crack tips. The strategy accelerates the convergence remarkably. Numerical experiments demonstrate the efficiency of the proposed algorithms applied to problems with several types of cracks.

Solution of large-scale 3D controlled-source electromagnetic modeling problem using efficient iterative solvers

With geophysical surveys evolving from traditional 2D to 3D models, the large volume of data adds challenges to inversion, especially when aiming to resolve complex 3D structures. An iterative forward solver for the controlled-source electromagnetic (CSEM) method requires less memory than that for a direct solver; however, it is not easy to iteratively solve an ill-conditioned linear system of equations arising from finite-element discretization of Maxwell’s equations. To solve this problem, we have developed efficient and robust iterative solvers for frequency- and time-domain CSEM modeling problems. For the frequency-domain problem, we first transform the linear system into its equivalent real-number format, and then introduce an optimal block-diagonal preconditioner. Because the condition number of the preconditioned linear equation system has an upper bound of [Formula: see text], we can achieve fast solution convergence when applying a flexible generalized minimum residual solver. Applying the block preconditioner further results in solving two smaller linear systems with the same coefficient matrix. For the time-domain problem, we first discretize the partial differential equation for the electric field in time using an unconditionally stable backward Euler scheme. We then solve the resulting linear equation system iteratively at each time step. After the spatial discretization in the frequency domain, or space-time discretization in the time domain, we exploit the conjugate-gradient solver with auxiliary-space preconditioners derived from the Hiptmair-Xu decomposition to solve these real linear systems. Finally, we check the efficiency and effectiveness of our iterative methods by simulating complex CSEM models. The most significant advantage of our approach is that the iterative solvers we adopt have almost the same accuracy and robustness as direct solvers but require much less memory, rendering them more suitable for large-scale 3D CSEM forward modeling and inversion.

Low-latitude reduction-to-the-pole and upward continuation between arbitrary surfaces based on the partial differential equation framework

SUMMARY Magnetic surveys conducted in complex conditions, such as low magnetic latitudes, uneven observation surfaces, or above high-susceptibility sources, pose significant challenges for obtaining stable solutions for reduction-to-the-pole (RTP) and upward-continuation processing on arbitrary surfaces. To tackle these challenges, in this study, we propose constructing an equivalent-susceptibility model based on the partial differential equation (PDE) framework in the space domain. A multilayer equivalent-susceptibility method was used for RTP and upward-continuation operations, thus allowing for application on undulating observation surfaces and strong self-demagnetization effect in a non-uniform mesh. A novel positivity constraint is introduced to improve the accuracy and efficiency of the inversion. We analysed the effect of the depth-weighting function in the inversion of equivalent susceptibility for RTP and upward-continuation reproduction. Iterative and direct solvers were utilized and compared in solving the large, sparse, non-symmetric and ill-conditioned system of linear equations produced by PDE-based equivalent-source construction. Two synthetic models were used to illustrate the efficiency and accuracy of the proposed method in processing both ground and airborne magnetic data. Aeromagnetic and ground data collected in Brazil at a low magnetic latitude region were used to validate the proposed method for processing RTP and upward-continuation operations on magnetic data sets with strong self-demagnetization.

Large-scale electromagnetic field analyses of coils wound with coated conductors using a current-vector-potential formulation with a thin-strip approximation

We developed a novel software for large-scale electromagnetic field analyses of coils wound with coated conductors based on current-vector-potential formulation with thin-strip approximation. Although this formulation was effective for obtaining the precise solutions of the electromagnetic field, the strong nonlinear property of superconducting materials frequently led to highly ill-conditioned linear systems of equations, which were difficult to solve efficiently. Moreover, the memory consumption and computation time required for the analyses rapidly increased with the size of the analysis due to dense matrix operations. In our software, the first difficulty was addressed by a novel preconditioning technique based on the algebraic multigrid method. Algebraic multigrid preconditioning enabled us to efficiently and stably solve the ill-conditioned linear systems of equations encountered in our analyses. It also improved the robustness of the analyses containing multifilament-coated conductors. As regards the second difficulty, the hierarchical matrices representation drastically reduced the memory consumption related to the dense matrices, as well as computation time. Meanwhile, our implementation of the hierarchical matrices representation was quite compatible with parallel computations on distributed memory computers. Finally, we presented some practical examples of large-scale analyses, which became possible by using the new software. For instance, the analysis of a cosine-theta dipole magnet whose number of degrees of freedom was more than 1.5 million was successfully completed in 78 h by 56 parallel processes and with a total memory consumption of 177 GB.

An improved iteration method for the numerical solution of groundwater flow in unsaturated soils

An improved Jacobi iterative method with two-side equilibration (TS-IJIM) is developed to model the process of unsaturated groundwater flow. The exponential model was adopted to linearize the nonlinear Richards equation (RE). Numerical discretization employing the finite difference method (FDM) produces a highly ill-conditioned system of linear equations. The TS-IJIM is proposed to reduce the matrix conditional number. Compared with two-side equilibration alone and the improved Jacobi iterative method alone, the proposed method effectively reduces the matrix conditional number, which is smallest when the resultant matrix has the same norm for each row. The TS-IJIM method is insensitive to the relaxation factor. Compared with the conjugate gradient method with a two-side equilibration algorithm, the TS-IJIM has better convergence in solving the highly ill-conditioned system. The conditional number of the ill-conditioned matrix can be reduced to 1.0 using the proposed TS-IJIM. Numerical results show that the proposed method can reduce the conditional number and improve convergence for solving unsaturated groundwater flow. Results also indicate that the proposed model is highly accurate and efficient in computing steady-state and transient 1D and 2D groundwater flow in unsaturated soils.

Parameter Identification in MHD Duct Flow Cauchy Problem using the DRBEM

In this study, the Hartmann number (Ha) is identified through the inverse formulation of Cauchy problem for magnetohydrodynamic (MHD) flow in a duct with insulated but no-slip or variably slipping walls. The solutions for direct and inverse problems are obtained using the dual reciprocity boundary element method (DRBEM), and Tikhonov regularization with L-curve method is used for solving ill-conditioned linear system of equations resulted from the inverse problem. The velocity and the induced magnetic field profiles are simulated from the direct solution depicting the well-known MHD characteristics. The Hartmann number is reconstructed from the inverse problem with an accuracy of 10−4 to 10−6. The DRBEM provides overdetermined information needed in the inverse MHD flow problem since it gives both the solution and its normal derivatives on the boundary points.

A dynamically deflated GMRES adjoint solver for aerodynamic shape optimization

An essential step in adjoint-based aerodynamic shape optimization is to obtain the Lagrange multiplier by solving a sparse linear adjoint system of equations based on the Jacobian matrices from the converged flow states. Such an approach has been applied widely within the aerospace community for the design of aircraft and other optimization problems for aerospace applications. However, the need to resolve the flow over complex geometries often requires highly stretched grids and gives rise to anisotropic flow fields which increase the stiffness of the discrete Jacobian needed for the solution of the adjoint system. When a GMRES algorithm preconditioned by an Incomplete LU factorization is used, this stiff linear system requires the use of a large number of Krylov subspace vectors and a high level of fill-in; both require an increase in the amount of memory. Deflated restarting, which distributes spectral eigenpairs, has proven to be an effective method to enhance the convergence rates when solving an ill-conditioned linear system of equations. In this paper, we propose an algorithm to dynamically determine the number of deflated vectors and to deactivate deflation to reduce memory utilization. We demonstrate its efficiency for a two-dimensional NACA 0012 airfoil, the three-dimensional Common Research Model (CRM) wing, and two three-dimensional complete aircraft configurations in viscous flow.

Block Diagonal Preconditioners for an Image De-blurring Problem with Fractional Total Variation

Even though using the total fractional-order variation (TFOV) model reduces or eliminates the stair casing effect in image de-blurring problems, but in the same time it generates a dense and ill-conditioned linear system of equations. These properties lead to slow convergence of any iterative method such as Krylov subspace methods. One treatment of the slowness property is to apply the preconditioning technique. In this research work, we propose two block diagonal preconditioners to construct an image from a given blurred image using circulant matrices. These matrices allow us to use the fast Fourier transform (FFT) and the convolution theorem. Using FFT and the convolution theorem reduces the cost of the computations from O(n2) into O(n log n) operations in each iteration and also reduces the storages. Our proposed preconditioners are of Murphy, Golub, and Wathen type. Numerical examples are given to illustrate the efficiency of our preconditioners.

Analysis of Basic Data of Unstable Problems of Coordinate Measurements of Geometrical Parameters of Products

The conditions for the appearance of instability have been analyzed in solving typical problems of coordinate measurements of the geometric parameters of products. The connection is shown between unstable problems of coordinate measurements and ill-conditioned systems of linear equations, and also with the parameters of the geometric elements that show signs of degeneracy. Formulas have been derived for estimating measurement error of the parameters of a circle, coaxiality of cylinders, and parallelism of planes. It has been established that the experimental results are consistent with the theoretical ones.

A novel Tikhonov regularization-based iterative method for structural damage identification of offshore platforms

The structural damage identification, which entails finding a best approximate solution to an over-determined system of linear equations, is actually one of the practical problems encountered in structural health monitoring. A major challenge to structural damage identification is the ill-conditioned problem caused by noise contamination and spatially incomplete measurements. A novel Tikhonov regularization iterative method (TRIM) is proposed to solve the ill-conditioned system of linear equations. This method iteratively reconstructs the regularization matrix by gradually updating the regularized solution. A merit is that the false-positive indicators of damage are greatly reduced; and as a result this approach would be able to detect smaller damages that could not be detected by  using traditional approaches. Another development embedded in TRIM is that a procedure, called singular value dichotomy, is developed to determine the regularization parameters. Several problems, such as the tremendous prior trials, due to the application of L-curve, are avoided. A numerical study is conducted on an offshore platform structure to demonstrate the effectiveness of the proposed method. The result shows that the new method outperforms the traditional Tikhonov regularization method in damage identification.

Efficiently solving total least squares with Tikhonov identical regularization

The Tikhonov identical regularized total least squares (TI) is to deal with the ill-conditioned system of linear equations where the data are contaminated by noise. A standard approach for (TI) is to reformulate it as a problem of finding a zero point of some decreasing concave non-smooth univariate function such that the classical bisection search and Dinkelbach’s method can be applied. In this paper, by exploring the hidden convexity of (TI), we reformulate it as a new problem of finding a zero point of a strictly decreasing, smooth and concave univariate function. This allows us to apply the classical Newton’s method to the reformulated problem, which converges globally to the unique root with an asymptotic quadratic convergence rate. Moreover, in every iteration of Newton’s method, no optimization subproblem such as the extended trust-region subproblem is needed to evaluate the new univariate function value as it has an explicit expression. Promising numerical results based on the new algorithm are reported.

On the Convergence of Iterative Methods and Pseudoinverse Approaches in Global Meshless Collocation

In the context of numerical approximation of partial differential equations, meshless methods are recent developments, which have been reported to be relatively straightforward, yet provide better convergence and accuracy as compared to the conventional mesh-based approaches, for some specific problems like stress–strain analysis and modeling in a deforming media. Among several proposed schemes, strong-form collocation schemes using radial basis functions are easy-to-program, require no mesh in the domain or at the boundary, avoid numerical integration, and have similar formulations for all dimensions. Some of the radial basiss functions like the Gaussian, multiquadric and inverse-multiquadric are infinitely smooth. The standard global meshless formulations based on such radial basis functions are often found to lead towards ill-conditioned systems of linear equations. For such formulations, we investigate the degree of ill-conditioning against degrees of freedom for different radial basis functions. Also, the convergence of four iterative methods and pseudoinverse approaches has been tested for the solution of such ill-conditioned systems. In order to compute the pseudoinverse, two different approaches, i.e., singular value decomposition and full rank Cholesky factorization have been used. A set of numerical experiments, performed here, demonstrate that pseudoinverse approach based on singular value decomposition performs better than the iterative methods. Although, pseudoinverse computation via full rank Cholesky factorization is faster, it is not recommended in global meshless collocation schemes due to poor convergence at relatively large degrees of freedom.

Numerical simulation of a one-phase steady flow towards a multistage fractured horizontal well

This paper presents a 3D mathematical model and its numerical implementation for a one-phase flow around a multistage fractured horizontal well with transverse fractures. The flows in the reservoir and in the fractures are governed by the Darcy’s law and were simulated separately. The problem was approximated by the finite volume method. The obtained systems of linear equations for the reservoir and the fractures were solved simultaneously, which allowed us to avoid the use of iterative process for solution adjustment both in the fractures and the reservoir. We investigated various techniques and proposed an optimal variant for solving ill-conditioned systems of linear equations appearing due to the grid approximation of the continuous mathematical flow model. We also studied how the fractures and well parameters affect the well productivity.

Multi-scale method for ill-conditioned linear systems

A multi-scale method is presented for solving ill-conditioned linear systems of equations. Numerical examples are provided to show the effectiveness and efficiency of the proposed method compared with traditional iterative methods.